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ABSTRACT 

We investigate the tidal interaction between a low-mass planet and a self- 
gravitating protoplanetary disk, by means of two-dimensional hydrodynamic sim- 
ulations. We first show that considering a planet freely migrating in a disk with- 
out self-gravity leads to a significant overestimate of the migration rate. The 
overestimate can reach a factor of two for a disk having three times the surface 
density of the minimum mass solar nebula. Unbiased drift rates may be obtained 
only by considering a planet and a disk orbiting within the same gravitational 
potential. In a second part, the disk self-gravity is taken into account. We con- 
firm that the disk gravity enhances the differential Lindblad torque with respect 
to the situation where neither the planet nor the disk feels the disk gravity. This 
enhancement only depends on the Toomre parameter at the planet location. It 
is typically one order of magnitude smaller than the spurious one induced by 
assuming a planet migrating in a disk without self-gravity. We confirm that the 
torque enhancement due to the disk gravity can be entirely accounted for by a 
shift of Lindblad resonances, and can be reproduced by the use of an anisotropic 
pressure tensor. We do not find any significant impact of the disk gravity on the 
corotation torque. 

Subject headings: accretion, accretion disks — hydrodynamics — methods: nu- 
merical — planetary systems: formation — planetary systems: protoplanetary 
disks 
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Introduction 



Since tlie discovery of the first exoplanet (IMayor &: Queloa Il995l ). theories of planet- 
disk interaction have receive d renewed attention. Using the analytic torque e xpression of 
Goldreich fc Tremaind (Il979l ) at Lindblad and corotation resonances, Ward (see IWardlll997l . 
and refs. therein) has elaborated a theory of planet-disk tidal interaction which shows that 
a planet embedded in a protoplanetary disk should experience an orbital decay toward the 
central object. For low-mass protoplanets, the timescale of this inward migration (usually 
known as type I planetary m igration) is much smaller than the disk lifetime, by typically one 
or two orders of magnitude (jWardI 119971 ). It puzzles current theories of planetary formation 
since it seems very unlikely that a giant planet can be built up before its protoplanetary core 
has reached the vicinity of the central star. 

Most of recent works dealing with planet-disk interac tions have therefore propose d mech- 
anisms that could slow down or stop type I migration. iMenou fc GoodmanI (120041 ) consid- 
ered realistic models of T Tauri a-disks instead of the customary power law models, and 
found that type I migration can be significantly slowed down at opacity transitions in the 
disk. iMasset et al.l (l2006bl ) showed that surface density jumps in the disk can trap low- 
mass protoplanets, thereby redu cing the type I migration rate to the disk's accretion rate. 
Paardekooper fc Mellemal (120061) found that the m i gratio n may even be reversed in disks of 
large opacity. More recently, iBaruteau fc Massetl (120081 ) have shown that, in a radiatively 
inefficient disk, there is an excess of corotation torque that scales with the initial entropy 
gradient at corotation. If the latter is sufficiently negative, the excess of corotation torque 
can be positive enough to reverse type I migration. 

A common challenge is in any case to yield precise estimates of the migration timescale. 
Nevertheless, a very common simplification of numerical algorithms consists in discarding the 
disk self-gravity. Apart from a considerable gain in computational cost, this is justified by 
the fact that protoplanetary disks have large Toomre parameters, so that the disk self-gravity 
should be unimportant. Even in disks that are not subject to the gravitational instability, 
neglecting the self-gravity may have important consequences on planetary migration, as we 
shall see. 

Thus far, a very limited number of works has taken t he di sk self-gravity into account 
in numerical simulations of planet-disk interactions. iBossI ( 120051 ) performed a large number 
of disk simulations in which the self-gravity induces giant planet formation by gravitational 
instability. His calculations are therefore short, running for a few dynamical times, and 
involve only very massive objects. The planets formed in these simulations excite a strongly 
non-linear res ponse of the disk, and a ny migration effects are probably marginal or negligible. 
Furthermore. iNelson fc Benz ( 2003al Jbl) included the disk self-gravity in their two-dimensional 
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simulations of planet-disk interactions. The authors find that the migration rate of a planet 
that does no t open a gap is slowed d own by at least a factor of two in a self-gravitating disk. 
Nonetheless, iPierens fc Hure J2005h (hereafter PH05) reported an analytical expression for 



the shifts of Lindblad resonances due to the disk gravity, and find that the disk gravity 
accelerates type I planetary migration. The apparent contradiction between these findings 
motivated our investigation. 

This work is the first part of a series of studies dedicated to the role of self-gravity on 
planetary migration. In the present paper, we focus on the impact of self-gravity on the 
migration of low-mass objects, that is on type I migration. This study will be extended 
beyond the linear regime in a future publication. 

The paper is organized as follows. The numerical setup used in our calculations is 
described in section |5J We study in section [3] the dependence of the differential Lindblad 
torque on the disk surface density, without and with disk self-gravity. We confirm in this 
section that the disk gravity accelerates type I migration, and check that this acceleration can 
be exclusively accounted for by a shift of Lindblad resonances. In section HI we show that the 
increase of the differential Lindblad torque due to the disk gravity can be reproduced with 
an anisotropic pressure tensor. We investigate in section [5] the impact of the disk self-gravity 
on the corotation torque. We sum up our results in section [61 



2. Numerical setup 

We study the impact of the disk self-gravity on the planet-disk tidal interaction by 
performing a large number of two-dimensional hydrodynamic simulations. Notwithstanding 
the need for a gravitational softening length, the two-dimensional restriction provides a direct 
comparison with the analytical findings of PH05 and enables us to achieve a wide exploration 
of the parameter space (mainly in terms of disk surface density, disk thickness and planet 
mass). 



2.1. Units 

As usual in numerical simulations of planet-disk interactions, we adopt the initial orbital 
radius Vp of the planet as the length unit, the mass of the central object as the mass unit 
and [GM^./rp^)''^/'-^ as the time unit, G being the gravitational constant (G = 1 in our unit 
system). We note Mp the planet mass and q the planet to primary mass ratio. 
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2.2. A Poisson equation solver for the code FARGO 



Our numerical simulations are performed with the code FARGO. It is a staggered mesh 
hydrocode that solves the Navier-Stokes and continuity equations on a po lar grid. It use s 
an upwind transport scheme with a harmonic, second-order slope limiter (Ivan Leerl 119771 ). 
Its particularity is to use a change of rotating fram e on each ring of the polar grid, which 
increases the timestep sig nificantly JMassetl EoOOal lbl) . thereby lowering the computational 
cost of a given calculation. 



2.2.1. Implementation 



We implemented a Poisson equation solver in FARGO as follows. Using the variables 
{u = logr, (p), where r and (p denote the polar coordinates, the potential V of the disk, as 
well as th e radial and azimuthal acc elerations Qr and g^p derived from it, involve convolution 
products (IBinney fc Tremaindll987l ). They can therefore be calculated at low-computational 
cost using Fast Fourier Transforms (FFTs), provided that a grid with a logarithmic radial 
spacing is used. Our Poisson equation solver calculates gr and g^ with FFTs. 

To avoid the well-known alias issue, the calculation of the FFTs is done on a grid 
whose radial zones number is twice that of the hydrodynamics grid, the additional cells 
being left empty of mass. Thus, the mass distribution of the hy drodynaniics ra esh can 
not interact tidally with its adjacent replications in Fourier space (ISellwoodI 119871 ) . and it 
remains isolated. Because of the the 27r— periodicity, such a precaution is not required in the 
azimuthal direction. 



Furthermore, a softening parameter e^g is adopted to avoid numerical divergences, the 
same way as the planet potential is smoothed. We point out that e^g must scale with r so 
that the expressions of gr and g^, smoothed over the softening length Ssg, involve indeed 
convolution products. The expressions of gr and g^p are given in Appendix |A1 

We finally present a test problem. For a two-dimensional disk with a uniform surface 
density S, gr reads 



gr{r)=4GJ: 



^(^^max) - K{V 

n 



EiUr 



'V. 



where K and E denote the complete elliptic integrals of the first and second kinds, respec- 
tively, where Wmin = r^^Jr and v^^^ = r/r^^^, r^m ('^max) denoting the disk inner (outer) 
edge (see PH05). We performed a self-gravitating calculation with S = 2 x 10~^, r^i 



and Tn 



2.5 Tp. The radial zones number is Nr 



min 0.4 Tp 

512, and we took a very small softening 
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length {esg{rp) is 100 times smaller than the grid radial spacing at r = Vp). Fig. [T] shows the 
agreement between the result of our calculation and the analytical expression of Eq. ([1]). The 
close-up displays gr around r = r^, for different softening length to mesh resolution ratios, 
e/Sr, at r = Vp. This shows the good convergence of our numerical calculation toward the 
analytical expectation when the softening length tends to zero. 

2.2.2. Numerical issues 

The implementation of the disk self-gravity addresses two issues. The first one concerns 
the convergence properties of our results. We performed preliminary runs to check the torque 
convergence, without and with self-gravity. The computational domain is covered with A^,. 
zones radially between r^i^ = 0.4rp and r^^x = 2.5rp, and zones azimuthally between 
ip = and (f = 2tt. For a comparative purpose, a logarithmic radial spacing is also used for 
the calculations without self-gravity. We adopted disk parameters and a planet mass that 
are representative of our study, namely a Q = 8 Toomre parameter at the planet location, 
and a Mp = 5 x 10~®M* planet mass. A complete description of our model parameters is 
deferred to section 12.31 We evaluate the torque obtained without self-gravity (Fnog) and 
with self-gravity (Ffsg) for several pairs {Nr,N^). The relative difference of these torques is 
displayed in Fig. [2^. We see in particular that the torque convergence is already achieved 
for Nr = 512 and = 1536, values that we adopted for all the calculations of this paper. 

Furthermore, since the softening length Egg varies from one ring to another, the FFT al- 
gorithm does not ensure an exact action-reaction reciprocity. Thus, the disk self-gravity may 
worsen the conservation of the total angular momentum (that of the system {gas-|-planet}). 
To investigate this issue, we performed calculations with a planet migrating in a disk without 
and with self-gravity. For these calculations only, the disk is inviscid, and reflecting bound- 
aries are adopted. As for the above convergence study, we adopted a Mp = 5 x 10~^M* 
planet mass, and a. Q = 8 Toomre parameter at the planet location. The value of £sg(^p) 
is the one used in our calculations hereafter (see section 12.31) . We display in Fig. the 
torques on the planet (Fpianct) and on the whole system (Fpianet+gas), for both calculations. If 
the code were perfectly conservative, the ratio Fpianct+gas/rpianct would cancel out, to within 
the machine precision. This ratio is typically ~ 0.5 % without self-gravity, and ~ 3 % with 
self-gravity. Although, as expected, the conservation of the total angular momentum is worse 
with self-gravity, it remains highly satisfactory. 
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Fig. 1. — Radial self-gravitating acceleration Qrir), in absolute value, for a uniform surface 
density field. The analytical expression of gr [see Eq. ([T])] is compared with the result of a 
self-gravitating calculation with a small softening length (see text). We point out that grir) 
is positive at the inner edge, then it becomes negative (here from r > 0.75). The close- 
up reveals the influence of the softening length on the agreement between the numerical 
calculation and the analytical expectation (see text). 
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2.3. Model parameters 

In the runs presented hereafter, the disk surface density S is initially axisymmetric 
with a power-law profile, S(r) = Ep(r/rp)~'^, where Ep is the surface density at the planet's 
orbital radius. The reference value of cr is 3/2. We therefore expect the corotation torque, 
which scales with the gradient of (the inyerse o f) the disk vortensity, to cancel out for a non 
self-gravitating disk JWardlll99ll : lMassetl[200ll ). 



The vertically integrated pressure p and S are connected by an isothermal equation 
of state, p = Scs^, where is the local isothermal sound speed. The disk aspect ratio 
is h{r) = H{r)/r = Cs{r)/rQK{f'), where H{r) is the disk scale height at radius r, and 
Qk denotes the Keplerian angular velocity. We take h uniform, ranging from h = 0.03 to 
h = 0.05, depending on the calculations. We use a uniform kinematic viscosity z/, which is 
10~^ in our unit system. 



The gravitational forces exerted on the disk include: 



The gravity of the central star. 

The gravity of an embedded planet, whose potential is a Plummer one with softening 
parameter e = 0.3H{rp). 

The disk self-gravity, whenever it is mentioned. The self-gravity softening length is 
chosen to scale with r, and to be equal to e at the planet's orbital radius, which yields 
£sg{r) = er/vp. Since h is taken uniform, H{r) scales with r, and e^ Jr) = 0.3H(r). 
We co mment that Bsg{rp) is very close to the recent prescription of iHure fc Pierens 



(120061 ) for the softening length of a fiat, axisymmetric self-gravitating disk. From now 



on, whenever we mention the softening length, we will refer to e. 



The disk's initial rotation profile Q{r) is slightly sub-Keplerian, the pressure gradient 
being accounted for in the centrifugal balance. When the disk self-gravity is taken into 
account, it reads 

nir) = (nlir) [I - (l + a)^] - '-^y\ (2) 

We comment that Qrir) is not necessarily a negative quantity. When it is so, the disk rotates 
slightly faster with self-gravity than without. In a two-dimensional truncated disk, Qr is 
positive at the inner edge and becomes negative at a distance from the inner edge that 
depends on a. We checked that, whatever the values of a used in this paper, gr is always 
negative in a radial range around the planet's orbital radius that is large enough to embrace 
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all Lindblad resonances (except the inner Lindblad resonance of m = 1 for a = 0, as can be 
inferred from Fig. [1]). 

As stated in section \2.2.2\ our calculations are performed on a grid with a logarithmic 
radial spacing, even when the disk self-gravity is not taken into account. The resolution 
is therefore the same in all our calculations. The computational domain is covered with 
Nr = 512 zones radially between r^i^ = OAvp and rmax = 2.5rp, and = 1536 zones 
azimuthally between ip = and (f = 27t. 



3. Dependence of the differential Lindblad torque on the disk surface density 

Our study is restricted to the linear regime, which enables us to compare the results of 
our calculations with analytical predictions. For this purpose, we consider a g = 5 x 10^^ 



planet to primary mass ratio. According to IMasset et al.l (j2006al ). for a two-dimensional 



calculation, the flow in the planet vicinity remains linear as long as 

tb < (3) 

where re = GMp/c^{rp) is the planet's Bondi radius and e is the softening length. Eq. ([3]) 
translates into q <ti qim, with = 0.3h^ in our units. For a h = 5 % disk aspect ratio. 



70 



gjin fa 4 X 10~ so that our planet mass is well inside the linear regime. For a. h = 3 
disk aspect ratio, gun ~ 8 x 10~^ and our planet mass approximately fulfills the linearity 
condition. Note that the linearity criterion given by Eq. ensures that the torque F exerted 
by the disk on the planet scales with g^. We focus in this section on the scaling of F with 
Sp, scaling that is expected, for a non self-gravitating disk, as long as the planet does not 



open a gap. The gap clearance criterion, recently revisited by iCrida et al.l (120061 ). reads in 
our unit system 

5.(|)-%50^,. ,4) 

The L.H.S. of Eq. (jl]) is ~ 100, hence we expect to check F oc Sp in our calculations without 
self-gravity. 

The runs presented hereafter lasted for 20 orbits, which was long enough to get sta- 
tionary values of the torque. For the calculations without self-gravity, the torque evaluation 
takes all the disk into account, it does not exclude the content of the planet's Hill sphere. 
We checked that excluding it or not makes no difference in the torque measurement. This is 
consistent with the fact that, for the planet mass considered here, we do not find any ma- 
terial trapped in libration around the planet, be it inside a circumplanetary disk (a fraction 
of the planet's Hill radius) or inside a Bondi sphere. 
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3.1. Case of a non self-gravitating disk 

We first tackle the case of a non self-gravitating disk. We measure the specific torque 
7 = r/g on the planet for six different values of Sp, ranging from Sp = 2 x 10~^ to Sp = 
2 X 10~^. This corresponds to varying the initial disk surface density at the planet's orbital 
radius from one to ten times the surface density of the minimum mass solar nebula (MMSN). 
Two situations are considered (see also Tabled]): 

• On the one hand, the planet does not feel the disk gravity: it is held on a fixed circular 
orbit, with a strictly Keplerian orbital velocity. In this case, referred to as the fixed 
case, both the planet and the disk feel the star gravity but do not feel the disk gravity. 
The disk non-Keplerianity is exclusively accounted for by the radial pressure gradient. 
This is the configuration th at has been contemplated in analytical torque estimates 



[see e.g. 



Tanaka et all 120021 ) 



• On the other hand, the planet feels the disk gravity. In other words, we let the planet 
evolve freely in the disk, so its angular velocity, which reads 

^pW = [^i'W-^7r.(rp)/rp]V2, (5) 

is slightly greater than Keplerian. In this case, which we call the free case, the planet 
feels the gravity of the star and of the disk while, as previously stated, the disk does 
not feel its own gravity. Contrary to the fixed case, the free case is not a self-consistent 
configuration since the planet and the disk do not orbit under the same gravitational 
potential. Nevertheless, this situation is of interest as it corresponds to the standard 
scheme of all simulations dealing with the planet-disk tidal interaction. 

Table 1. Planet's angular velocity ^lp{rp) and disk's rotation profile Q{r) for a non 

self-gravitating disk 



fixed case free case 



2 ,„ ^ 9r(rp)\l/2 



Qpirp) ^K{rp) (njf{rp) 

n{r) nKir)[l-{l + a)h^]^^^ C^W [l - (1 + 

In both cases, the initial planet's angular velocity is strictly Ke- 
plerian 

For all the runs presented here, gr(fp) < so that Qp{rp) is 
slightly greater in the free case than in the fixed case 
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We show in Fig. [3] the specific torques (in absolute value) obtained with the fixed and 
free cases, for a. h = 0.05 disk aspect ratio. In the fixed situation, there is an excellent 
agreement with the expectation 7 oc Sp, and, not surpr isingly, the to r ques a re bounded by 



the two- and three-dimensional analytical estimates of iTanaka et al.l (|2002| ). Nonetheless, 
the free case reveals two unexpected results. For a given surface density, the absolute value 
of the torque is larger than expected from the fixed case. Moreover, it increases faster than 
linearly with the disk surface density. 

The two latter results can be explained with the relative positions of the Lindblad 
Resonances (hereafter LR) in the fixed and free cases. We display in Fig. the locations 
Hlr ('"olr) of an Inner (Outer) LR, when the planet is on a fixed orbit. They are given by 
Hlr = ^^^(^ilr) and tqlr = ^"^(^olr), with Q{r) the disk's rotation profile (solid curve), 
and fiiLR (^^olr) the frequency of the ILR (OLR), simply deduced from the planet frequency 

When the planet is on a free orbit (Fig. |Dd), its frequency is slightly larger than in the 
fixed case. Thus, the frequencies of the LR are also larger in the free case, which induces a 
spurious inward shift of all the resonances. The OLR get closer to the orbit, which increases 
the (negative) outer Lindblad torque. The ILR are shifted away from the orbit, which 
reduces the (positive) inner Lindblad torque. Thus, the (negative) differential Lindblad 
torque is artificially larger in the free case. 

The inward shift of the LR, which we denote by 6R, has been evaluated analytically 
by PH05. A simple estimate can be obtained as follows. We denote by i?* the nominal 
position of the resonances without disk gravity. We assume that the disk's rotation profile 
is strictly Keplerian. The shift 6R being induced by the increase of the planet frequency, we 
have SR/R^ = —26flp{rp)/3flK{^p), where 6flp{rp) is the difference of the planet frequencies 
between the free and fixed cases. Using Eq. and a first-order expansion, we are left with 

R* 3rpr2^ (^p) 

A more accurate expression for 5R/R^. is given by PII05 [see their equation (7c)]. Eq. ([6]) 
shows that the shift of the LR scales with Qrifp), hence with Sp. This explains why the 
torque in the free case increases faster than linearly with the disk surface density. The 
relative shift of the resonances 5R/R^ typically amounts from —3 x 10~^ to —3 x 10"'^ for 
our range value of surface densities, corresponding however to a torque relative discrepancy 
between ~ 12 % and ~ 120 % (see Fig. [3]). 

We are primarily interested in a quantitative comparison of the torques in the fixed and 
free cases. Nonetheless, since the shift of the LR scales with gr{rp), it depends on the mass 
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distribution of the whole disk. Thus, the torque discrepancy between the fixed and free cases 
also depends on gr{rp), hence on Sp, a, rmin and Tmax- In particular, we point out that if the 
planet is close enough to the disk's inner edge, then gr{rp) can be positive (see Fig. [H for 
cr = 0). This shifts all the LR outward (instead of inward) and reduces the torque. We have 
checked this prediction with an appropriate calculation (not presented here). 

In our study, only Sp is a free parameter. The index of the unperturbed surface density 
profile, a, is fixed indeed to 3/2, as explained in section 12. 3[ Our values of Vjam and r^^.^ 
are those customarily u sed in numerical simulations of planet-disk interactions (see e.g. 



de Val-Borro et al.l 120061 ). Thus, a useful quantitative comparison of the torques between 
the free and fixed cases can be provided just by varying Ep. In particular, one may think 
the torque discrepancy to be significant only for high values of Sp. Nevertheless, such a 
discrepancy depends both on the surface density Sp and on the disk aspect ratio h. As 
explained in Appendix [Bl we expect the relative difference of the torques between the free 
and fixed situations to scale with (Qh)"^, where Q is the Toomre parameter at the planet's 
orbital radius. 



with K the horizontal epicyclic frequency, defined as k = [2Qr~'^ d{r'^Q) / drf'^'^ , and mo = 
TTVpTip/M^. Eq. ([7]) can be recast as Q = /i/7rSp in our units. 

To study the impact of h on previous results, we performed another set of calculations 
with h = 0.03. Fig. [5] confirms that the relative difference of the torques scales with the 
inverse of Qh. It yields an estimate of the error done on the torque evaluation when involving 
the strongly biased free situation rather than the self-consistent fixed situation. For instance, 
for a h = 3 % disk aspect ratio, the free situation can overestimate the torque by as much 
as a factor two in a disk that has only ~ 3 times the disk surface density of the MMSN. 
Moreover, the torque relative difference is less than 20 % as long as Qh > 2.5, hence as 
long as the Toomre parameter at the planet location is approximately greater than 50 if 
h = 0.05, or 80 if = 0.03. Remember that these estimates depend on the precise value of 
gr{rp), hence on the mass distribution of the whole disk. They are provided with fixed, but 
customarily used values of a, rmin and rmax- 

To avoid the above torque discrepancy, one must ensure that the planet and the disk 
feel the same gravitational potential. The workaround depends on whether the disk is self- 
gravitating or not, and whether the planet freely migrates in the disk or not: 



1. The disk is not self-gravitating. The planet's angular velocity should therefore be 
strictly Keplerian: 
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(a) The planet evolves freely in the disk. Thus, its angular velocity, given by Eq. ([5]), 
is slightly greater than Keplerian. A workaround could be to subtract the ax- 
isymmetric component of the disk surface density to the surface density before 
calculating the force exerted on the planet by the disk. This would cancel out 
Qrifp), and the planet's angular velocity would remain strictly Keplerian. 

(b) The planet is held on a fixed circular orbit, with necessarily a Keplerian angular 
velocity. This is a self-consistent situation. 

2. The disk is self-gravitating. The planet's angular velocity should therefore be given by 



(a) The planet evolves freely in the disk. This is a self-consistent situation. 

(b) The planet is held on a fixed circular orbit. This situation is self-consistent only 
if the planet's fixed angular velocity is given by Eq. ([5]). 

From now on, whenever calculations without disk gravity are mentioned, they refer to 
the fixed situation. We mention them as nog calculations. 



We study how the results of section 13.11 differ when the disk gravity is felt both by the 
planet and the disk. The planet is still held on a fixed circular orbit at r = rp, its angular 
velocity is given by Eq. ([5]). As in the situation without disk gravity, the planet's initial 
velocity is that of a fluid element that would not be subject to the radial pressure gradient 
(see Table E]). 

Taking the disk self-gravity into account induces two shifts of Lindblad resonances 
(PH05): (i) a shift arising from the axisymmetric component of the disk self-gravity, and 
(ii) a shift stemming from the non-axisymmetric component of the disk self-gravity. We 
therefore performed two series of calculations: 

Table 2. Planet's angular velocity fip(rp) and disk's rotation profile fl{r), without and 

with disk gravity 



Eq. (ISD: 



3.2. Case of a self-gravitating disk 



Without disk gravity 



With disk gravity 




) 



1/2 
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1. Calculations that involve only the axisymmetric part of the disk self-gravity. They 
are mentioned as axisymmetric self-gravitating calculations (asg calculations). Their 
computational cost is the same as that of a calculation without disk gravity since only 
one-dimensional FFTs are performed. The results of these calculations are presented 
in section [3. 2. 1[ 

2. Fully self-gravitating calculations (fsg calculations), which are more computationally 
expensive as they involve two-dimensional FFTs. Their results are presented in sec- 
tion [3D 

3.2.1. Axisymmetric self-gravitating calculations 

We display in Fig. [6] the torques obtained with the nog, asg and fsg calculations, when 
varying Sp. We will comment the results of the fsg calculations in section [3.2.21 The torques 
obtained in the asg situation, which we denote by 7asg, are hardly distinguishable from the 
torques without disk gravity, mentioned as 7nog- A straightforward consequence is that 
7asg scales with Sp with a good level of accuracy. We point out however that the torque 
difference |7asg| — l7nog| is slightly negative and decreases with Sp (not displayed here). The 
relative difference | |7asg| — |7nog| |/|7nog| varies from ~ 0.2 % for Ep = 2 x 10~^, to ~ 2 % for 
Sp = 2 X 10-3. 

The interpretation of these results is as follows. In the asg situation, the positions of 
the LR related to the Fourier component with wavenumber m are the roots of equation (see 
PH05 and references therein) 

^asg(r) = /t'(r) - m^[Vl{r) - Vl^f + m^cl{r)/r^ = 0, (8) 

where, contrary to the nog situation, VL{r) and Vtp depend on gr (see Table [2]). As in 
section 13.11 the increase of the planet frequency implies an inward shift of the LR, which 
increases the differential Lindblad torque (see Fig. [Dd). Furthermore, as pointed out in 
Fig. [7^1, the increase of the disk frequency causes an outward shift of all LR, which reduces 
the differential Lindblad torque. Accounting for the axisymmetric component of the disk 
gravity therefore leads to two shifts of the resonances, acting in opposite ways. Fig. [Tb 
shows that both shifts do not compensate exactly: the LR are slightljfl moved away from 
corotation with respect to their nominal position without disk gravity. This is in qualitative 



^To improve the legibility of Figs. [3| and [3 the disk's rotation profile with self-gravity is depicted with a 
value of Ep that is 25 times greater than the maximal value of our set of calculations. 
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agreement with PH05, who found a resulting shift which is negative for inward resonances, 
and positive for outward resonances (see their 6Ri + 6R3 expression). The sign of the shift 
results from the fact that the disk's rotation profile decreases more slowly with self-gravity 
than withoul§|, and explains why |7asg| — l7nog| is negative. The absolute value of this shift 
increases with Sp, which entails that ||7asg| — |7 nog 1 1 increases with Sp. 



3.2.2. Fully self- gravitating calculations 

We now come to the results of the fsg calculations depicted in Fig. El The torques 
obtained with the fsg calculations, denoted by 7fsg, are larger than 7asg and 7nog- Moreover, 
7fag;| grows faste r thari linearly with the disk surface density, a result already mentioned by 



Tanigawa fc LinI (120051 ). 



These results can be understood again in terms of shifts of the LR. Besides the shift due 
to the slight increase of the planet and of the disk frequency, the fsg situation triggers another 
shift stemming from the additional non-axisymmetric term — 27rGEm/r in the dispersion 
relation of density waves (in the WKB approximation, see PH05). The positions of the LR 
associated with wavenumber m are this time the roots of equation 

Dfsg(r) = Dasg(r) - 27rGS(r)m/r = 0, (9) 

where Dasg is given by Eq. (jS]). PH05 showed that: 

• This non-axisymmetric contribution moves inner and outer LR toward the orbit, with 
respect to their location in the asg situation. This explains why |7fsg| > |7asg|, and 
implies that the torque variations at inner and outer resonances have opposite signs. 

• The shift induced by the non-axisymmetric part of the disk self-gravity dominates that 
of its axisymmetric component. Therefore, it approximately accounts for the total shift 
due to the disk gravity, and explains why |7fsg| > |7nog| ~ l7asg|- 

• This shift increases with Sp, so that |7fsg| increases faster than linearly with Hp. 

Our results of calculations are in qualitative agreement with the analytical work of PH05. 
Before coming to a quantitative comparison in section 13. 3. 2^ we focus on the relative differ- 
ence of the torques between the fsg and nog situations. From previous results, we assume 



^We comment that this statement is not straightforward since it involves both the sign and the variations 
of function gr] here again we checked that this statement is vahd in a radial range around the orbit that is 
large enough to concern all LR. 
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that the only shift of the LR is due to the non-axisymmetric part of the disk gravity. In- 
terestingly, this shift does not feature Qr, so it does not depend on the mass distribution 
of the whole disk. It only depends on the surface density at the planet location. Since the 
torque variations at inner and outer resonances are of opposite sign, we expect the relative 
difference of the torques to scale with Q~^, for high to moderate values of Q. This is shown 
in Appendix O It differs from the {Qh)~^ scaling obtained in Fig. [5l where the torque 
variations at inner and outer resonances were of identical sign. 

In Fig. [HI we plot this relative difference as a function of Q for previous results and 
for another series of runs performed with a, h = 0.03 disk aspect ratio. The departure 
from the expected scaling occurs for (5^6. The behavior at low Q will be tackled in 
section 13.3.21 Fig. [S] yields a useful estimate of the torque increase due to the disk gravity, 
or, differently stated, of the torque underestimate if one discards the disk gravity. As such 
estimate only depends on the Toomre parameter at the planet location, whatever the global 
mass distribution of the disk. The torques' relative difference is typically one order of 
magnitude smaller than in the situation of a planet freely migrating in a non self-gravitating 
disk (Fig. [5]). It amounts typically to 10 % for Q «i 10. For Q > 50, accounting for the disk 
gravity or not has no significant impact on the torque measurement. 

Our results confirm that the disk gravity acce l erates t ype I migration. This might sound 



contradictory with the results of iNelson &: Bend (l2003a| ]bl). who found that the disk self- 



gravity slows down migration for a planet that does not open a gap. The authors compared 
however the results of their self-gravitating calculations (where both the planet and the 
disk feel the disk gravity) to those obtained with the misleading situation of a planet freely 
migrating in a disk without self-gravity. As shown by Fig. [9l or as can be inferred from Figs. [3] 
and [6l comparing both situations would lead us to the same conclusion. There is therefore 
no contradiction between their findings and ours. From now on, we do not distinguish the 
gravity and self-gravity designations, since the planet and the disk orbit within the same 
potential in our calculations. Whenever calculations with disk gravity are mentioned, they 
refer to the fsg situation. 



3.3. Comparison with analytical results 

3.3.1. An analytical estimate 

We propose in this section a simple analytical estimate of the relative difference of 
the torques between the fsg and nog situations. This estimate concerns high to moderate 
values of the Toomre parameter at the planet location. We assume that the only shift of 
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the LR in the fsg situation arises from its non-axisymmetric contribution. This comes to 
approximating the nog and asg situations, which is a reasonable assumption from Fig. [61 
Furthermore, since this shift has same order of magnitude at inner and outer LR (PH05), we 
focus on the one-sided Lindblad torque and use a local shearing sheet approximation. We 
set up local Cartesian coordinates (x, y) with origin at the planet position, the x and y-axis 
pointing toward the radial and azimuthal directions. Our x-coordinate is taken normalized 
as 

(10) 



r — Tr. 



r — r.n 



X 



As is usually done in the shearing shee t framework, we disc ard the radial dependence of the 
disk surface density and scale height ( iNarayan et al.l 119871 ). In a non-gravitating disk, the 
LR associated with wavenumber m are therefore located at 



X 



nog 



2 y/T+e 



where ^ = mh, e = +1 for outer resonances, e = — 1 for inner resonances. In the fsg situation, 
LR are located at Xnog + Sxfsg, where the shift 6xisg is evaluated by -Dfsg(xnog + ^a^fsg) = 0. 



Using Eqs. 



and ffTT]) . a first-order expansion yields 

2 1 



SXf: 



3eQ ^/TTe 



(12) 



We comment that the equation (7b) of PH05 reduces to our Eq. (fT2|) for a surface density 
profile decreasing as r~^/^. 

In the linear regime, the one-sided Lindblad torque F amounts to a summation over 
m of the Fourier components F^^. In the shearing sheet approximation, since all quantities 
depend on m through ^, the summation over m is approximated as an integral over C,, 



-1 



T{x = xl,0 d^, 



(13) 



where xl denotes the positions of t he LR, T is the m}^ Fourier component of the one-sided 
Lindblad torque, given by (see e.g. IWardlll997l ) 



Tix,0=K- 



(14) 



with K a constant. We assume that Eq. IHM can be used whatever the disk is self-gravitating 
or not ( Goldreich fc TremainelE979l ). The forcing function \& in Eq. 0141) is approximated in 
a standard way as a function of the Bessel Kq and Ki functions. 



v^(x, = i^i(|x|0 + 2 VI + e^o(|a;|0 



(15) 
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We furth ermore approximate ^ (x, as (|a;|.^)~^, to within a numerical factor of the order 
of unity (lAbramowitz fc Stegunlll972l ). This approximation is vahd when ^ ^ 1, hence for 
low m- values. 

With a first-order expansion in Q~^, the difference of the one-sided Lindblad torques 
between the fsg and nog situations reads 



h 



Combining Eqs. (fTTj) to (fT6l) . we are left with 

Tfsg -^nog 



where 



3 X 



21 



e 



(i+a^/'(i+4e 



(1 + ^^)3/2(1+4^2 

2^3 - log (7 + 4^3) 



V3-n/3 



1.21. 



(16) 



(17) 



Not surprisingly, the relative difference of the one-sided Lindblad torques scales with the 
inverse of Q. This is the same scaling as for the relative difference of the differential Lindblad 
torques, assuming high to moderate values of Q (see Appendix O and Fig. [8]). Note that, 
unlike the analysis of PH05, the present analysis, which is restricted to the shearing-sheet 
framework, enables one to exhibit the scaling given by Eq. ( JT7j) . 



3.3.2. Results 

We come to a quantitative comparison of our results of calculations with our analytical 
estimate, given by Eq. (fT7|) . and the analytical results of PH05, who estimated the depen- 
dence of the differential Lindblad torque on the disk mass, for a fully self-gravitating disk 
(see their figure 4b). Another series of fsg calculations was performed with disk parameters 
similar to those of PH05, namely a. h = 5 % disk aspect ratio, a planet mass corresponding 
to the linear regime (its value is precised hereafter). We vary the disk surface density at 
the planet's orbital radius from Sp = 4 x 10~^ to Sp = 10^^, This corresponds to varying 
Q from 40 to 1.6. The runs lasted for 10 planet's orbital periods, which was long enough 
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to get stationary torques for the largest values of Q, but short enough to avoid a signifi- 
cant growth of non-axi symmetric per turbations for the lowest values of Q, probably due to 



SWING amplification flToomrelll964f ) 



As we aim at comparing the results of two-dimensional calculations with analytical 
expectations (for which there is no softening parameter), we investigated how much our 
results of calculations depend on the softening length. For this purpose, the calculations 
were performed with three values of e: 0.1H{rp), 0.3H{rp) and 0.6H{rp). The planet mass 
is Mp = 4.4 X lO-^M, for e = 0.3H{rp) and e = 0.6H{rp), whereas Mp = IQ-^M, for 
e = 0.1H[rp). This choice ensures that the Bondi radius to softening length ratio does not 
exceed ~ 10 % for each value of e. 

Each calculation was performed with and without disk gravity, so as to compute the 
relative difference of the torques between both situations. The reason why we compute this 
relative difference is that it does not depend on the details of the torque normalization, 
be it for the numerical or the analytical results. Nonetheless, PH05 only calculated the 
normalized torque in the fsg situation as a function of the disk mass. We then evaluated 
their normalized torque without disk gravity by extrapolating their torque with disk gravity 
in the limit where the disk mass tends to zero. 

Fig. [To] displays the relative difference of the torques between the fsg situation (7fsg) 
and the nog situation (7nog); obtained with our calculations, the analytical expectation of 
PH05 and our analytical estimate. This relative difference grows faster than linearly with 
Sp, although a linear approximation is valid at low surface density, as already stated in 
section 13.2.21 Our linear estimate is in agreement with the results of calculations with 
e = 0.6H{rp) up to Q ~ 3, where it leads to a torque enhancement that is typically half 
the one estimated by PH05. Furthermore, our results of calculations depend much on e, 
more especially at high Ep. For a given value of Ep, the relative difference of the torques 
decreases as e increases. Differently stated, increasing the softening length reduces 7fsg more 
significantly than 7nog- 

We finally comment that our results of calculations with e = O.lif(rp), which matches 
the mesh resolution in the planet vicinity, are in quite good agreement with the analytical 
prediction of PH05. Surprisingly, the relative differences obtained with our calculations are 
greater than their analytical expectation. We checked that doubling the mesh resolution 
in each direction does not alter the relative differences measured with our calculations, 
as already pointed out in section 12.2.21 (see Fig. [2^). We show in Appendix [D] that this 
result can be explained by the failure of the WKB approximation for low values of the 
azimuthal wavenumber. The relative difference between the results of our calculations and 
the predictions of PH05 is ~ 15 % for Q ~ 8, and does not exceed ~ 30 % for Q < 2. 
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This satisfactory agreement confirms that the impact of the disk gravity on the differential 
Lindblad torque may be exclusively accounted for by a shift of Lindblad resonances. 

4. Modeling the non-axisymmetric contribution of the disk self-gravity with 

an anisotropic pressure tensor 

In section [3|, we investigated the impact of the disk gravity on the differential Lindblad 
torque for low-mass planets. The torque of an asg calculation (where only the axisymmetric 
component of the disk self-gravity is taken into account) is close to that of a nog calculation 
(without disk gravity). However, a fsg calculation (which furthermore involves the non- 
axisymmetric contribution of the self-gravity) displays a significant increase of the torque, 
which can be exclusively accounted for by a shift of the LR. 

We propose in this section to model this torque enhancement for low-mass planets. Our 
model aims at calculating only the axisymmetric part of the disk self-gravity, and applying 
an additional shift of the LR that mimics the one of its non-axisymmetric part. Altering 
the location of the LR comes to modifying the dispersion relation of the density waves. The 
dispersion relations of the asg and fsg cases differ only from the — 27rG'Sm/r term [in the 
WKB approximation, see Eqs. ([H]) and (jH])]. There is however no straightforward way to 
add an extra term proportional to m in the dispersion relation -Dasg of the asg situation. 
We propose to multiply the m'^cl/r^ term of Z^asg by a constant, positive factor 1 — a, with 
a > to ensure that LR are shifted toward the orbit. This can be achieved by multiplying 
the azimuthal pressure gradient d^P by 1 — a in the Navier-Stokes equation or, differently 
stated, by assuming an anisotropic pressure tensor, for which the pressure in the azimuthal 
direction reads Pip = (1 — C()Pr, where Pr, the pressure in the radial direction, is given 
by Pr = Sc^. We call a the anisotropy coefficient. When an asg calculation includes the 
anisotropic pressure model, it is mentioned as an asg-|-ap calculation. We comment that the 
rotational equilibrium of the disk, which involves the radial pressure gradient, is not altered 
by this model. 

We now explain how to take the adequate value for the anisotropy coefficient. As in 
section |3l we assume an initial surface density profile scaling with r~^/^, inducing a negligibl^ 
vortensity gradient, hence a negligible corotation torque. Thus, the torques obtained with 
our calculations only include the differential Lindblad torque. We denote by Ffgg, Fasg and 
Fasg+ap the differential Lindblad torques obtained with the fsg, asg and asg+ap calculations. 



■^With a uniform disk aspect ratio, the vortensity gradient vanishes for a non self-gravitating disk while 
it is negligible, but does not cancel out, for a self-gravitating disk. 
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Our model aims at imposing that 

Tasg+ap -Tasg -^fsg Tasg- (-^'^) 

A first-order expansion of tlie L.H.S. of Eq. ( |T9l) witli a, and of its R.H.S. witli Q^^ leads to 

a = PQ-\ (20) 

where 

(drasg+ap/C'a)a=0 

The parameter j3 depends only on the softening length to disk scale height ratio rj = e/H. 
We calculated it for rj = 0.1, 0.3 and 0.6 for small, fixed values of a and Q^^, which we 
denote with a zero subscript. For each value of r], we performed an asg, an asg+ap and a fsg 
calculation with q = 10~^ and h = 5 %, corresponding to a Bondi radius to softening length 
ratio of ~ 2.7 %. Furthermore, we adopted Sp = 5 x 10""^, yielding Q^^ ~ 0.03. The asg+ap 
calculation had = 0.01. Using Eq. (12T|) . the parameter (3 was therefore calculated by 

P = aoQo /''''^';' ■ (22) 

^ asg+ap ^ asg 

We display in Table [3] the values of P for r] = 0.1, 0.3 and 0.6. We note that our anisotropic 
pressure model should be applied only when Q > P to satisfy the constrain 1 — « > 0. This 
is not a stringent constrain since (3 < 1 for these values of rj. 

We comment that the value of m for which the resonance shifts induced by the self 
gravity and by the anisotropic pressure are equal is beyond the torque cut-off. Several 
reasons may conspire for that: 

• For a given shift, the relative torque variation is larger for resonances that lie closer to 
the orbit, which gives more weight to high— m component. 

• The shifts estimated by a WKB analysis may dramatically differ from the real shifts 
(see Appendix [D]) , especially at low— m, where significant torque is exerted. 

• The torque expression for an anisotropic pressure h as not been worked out in the liter- 



ature, and may differ from the standard expression (jWardlll997l ). with the consequence 



that equal shifts will not yield equal torque variations. 
Table 3. Calculation of the anisotropy coefficient: values of (3 for different values of rj 



ri = e/H 0.1 


0.3 


0.6 


/3 0.32(4) 


0.61(4) 


0.94(1) 
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4.1. Validity of the anisotropic pressure model 

We first test the validity of our model by performing a series of calculations with Q 
ranging from 1.5 to 8. From Eq. ([7]), Q can be set by varying either /z or Ep. Varying h 
however alters the ratio r^/e, which controls the flow linearity in the planet vicinity. We 
therefore fixed h = 0.05 and varied Sp. The planet to primary mass ratio is g = 10~^, the 
softening length is e = 0.3H{rp). For each value of Q, we performed a fsg, an asg, and an 
asg+ap calculation, for which the anisotropy coefficient is a = P/Q, with P = 0.614 (see 
Table [3]). The results are displayed in Fig. [TTJ As expected from the first-order expansion 
in used to derive Eq. (12U]) . the difference between the torques of the fsg and asg+ap 
calculations increases when decreasing Q. The relative difference is ~ 0.4 % for Q = 8, 
~ 10 % for Q = 2.5, and reaches ~ 25 % for Q = 1.5. The anisotropic pressure model 
therefore reproduces the torque of a fsg calculation with a good level of accuracy up to 
Q~4. 

The robustness of our model is furthermore tested against the onset of non-linearities, 
by varying the planet to primary mass ratio q. The Toomre parameter at the planet location 
is fixed at Q = 8. A series of asg, asg-fap and fsg calculations was performed with q ranging 
from 10^® to 7 X 10~^, hence with tb/b ranging from ~ 2.7 % to ~ 18.7 %. Fig. [T2] displays 
the specific torque as a function of q for each calculation. The torques obtained with the 
fsg and asg+ap agree with a good level of accuracy. Their relative difference, shown in the 
close-up, increases almost linearly from ~ 0.4 % to ~ 4 %, due to the onset of non-linearities. 

These results indicate that the anisotropic pressure model succeeds in reproducing the 
total torque obtained with a fully self-gravitating disk, as far as a low- mass planet, a high to 
moderate Toomre parameter, and a surface density profile scaling with r~^/^ are considered. 
With these limitations, these results present another confirmation that the impact of the disk 
gravity on the differential Lindblad torque can be entirely accounted for by a shift of the 
LR. We suggest that in the restricted cases mentioned above, the anisotropic pressure model 
could be used as a low-computational cost method to model the contribution of the disk 
gravity. We finally comment that, not surprisingly, these results do not differ if the planet 
freely migrates in the disk, which we checked with long-term fsg and asg+ap calculations 
(not presented here). 

5. Corotation torque issues 

Hitherto, we have considered an initial surface density profile scaling with r~^/^, induc- 
ing a negligible vortensity gradient, hence a negligible corotation torque. This assumption 
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ensured that the torques derived from our calculations accounted only for the differential 
Lindblad torque. It enabled a direct comparison with analytical expectations focusing on the 
differential Lindblad torque. We release this assumption and evaluate the impact of the disk 
self-gravity on the corotation torque Fc, in the linear regi me. For a disk witho ut self-gravity, 





(Ward 


1991. 


1992; 


Masset 


2001) 



dln(S/E) 
dXnr 



(23) 



where Xs is the half-width of the horseshoe region, Tc denotes the corotation radius, and 
B = {2r)^^ dijr'^Vt) / dr is half the vertical component of the flow vorticity. We denote by 
rc,asg5 Tcasg+ap, and Fc,fsg the corotation torques in the asg, asg+ap, and fsg situations. 
The same quantities without the C subscript refer to the total torque in the corresponding 
situation. 



We performed the same set of asg, asg+ap and fsg calculations as in section 14. H but 
with a flat initial surface density profile (we vary the planet to primary mass ratio g, for 
Q = 8). An additional nog calculation was also performed for g = 5 x 10~®. The results of 
these calculations are displayed in Fig. [121 The torques of the nog and asg calculations are 
hardly distinguishable, their relative difference being ~ 2 %, similarly as in section 13.2. ![ 
where a = 1.5. This difference should therefore be attributed to the differential Lindblad 
torque rather than to the corotation torque. It confirms that the corotation torque is not 
altered by the axisymmetric component of the disk gravity. 

Furthermore, the torques of the fsg runs are significantly larger than those of the asg+ap 
runs. Their relative difference varies from ~ 11 % to ~ 17 %. We do not expect this difference 
to arise from the differential Lindblad torque, despite the change of a. The differential 
Lindblad torques should therefore differ from ~ 0.4 % to ~ 4 %, as for a = 1.5 (close-up of 
Fig. [12]). This reveals that the fsg situation, or the asg+ap situation, or both, boosts the 
(positive) corotation torque. 



We expect in fact the asg+ap situation to enhance the corotation torque. iMasset et al. 



(l2006al ) have estimated Xg for a disk without self-gravity, in the linear regime. Their estimate 
reads Xg ~ l.lQrp^/qJh. In the limit where the planet mass vanishes, a fluid element on a 
horseshoe separatrix has a circular trajectory and is only sensitive to the azimuthal gradient 
of the disk pressure. The above estimate of Xg therefore holds for an asg+ap calculation 
if one substitutes h with a/1 — a h, which we checked by a streamline analysis. Thus, we 
expect the anisotropic pressure model to slightly increase the half-width of the horseshoe 
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zone, thereby increasing the corotation torque as 

with Fcasg given by Eq. (!23|) . and a = [3/Q. 

To investigate whether the fsg situation also increases the corotation torque, we evaluate 
the quantity (Fcfsg — rc,asg)/rc,asg5 which can be recast as 

Tcfsg TCjasg TCifsg r^^asg+ap TCjasg+ap TCjasg (0^\ 

f ~ f f ■ ^ ' 

-•- C,asg J- C,asg J- C,asg 

Using Eq. fl2i|) . the second R.H.S. of Eq. fl25|) reads a /{I — a), and is ~ 8.4 %. Moreover, 
for the sake of simplicity, we neglect the relative change of the differential Lindblad torques. 
This assumption is grounded for the smallest planet masses that we consider, for which, as 
stated above, this change does not exceed ~ 1 %. The first R.H.S. of Eq. therefore reads 
(F fag ~ rasg+a.p)/rc,a,sg . The quantity Fcasg can be connected with Fasg, using the estimate 
of 



Tanaka et al.l (120021 ) for a flat surface density profile. This connection is motivated by 
the fact that both the differential Lindblad torque, and the corotation torque are almost 
identical in the nog and asg situations. This leads to Fcasg ~ — l-56Fasg. Eq. fl25|) finally 



reads 



Fc,fsg Fc^asg -^fsg Fasg+ap 01 



"C,asg 1.56 Fasg I- a 



+ q • (26) 



This ratio is displayed in the close-up of Fig. [131 It shows that the fsg situation slightly 
enhances the corotation torque, but this enhancement does not exceed ~ 4.5 % for the 
highest planet mass that we consider. For the smallest planet masses, it is negligible with 
respect to the increase of the corotation torque triggered by the asg+ap situation. Thus, the 
large difference between the torques of the asg+ap and fsg calculations can be exclusively 
accounted for by the boost of the corotation torque in the asg+ap situation. 

The slight increase of the corotation torque in the fsg calculations should be compared 
to that of the differential Lindblad torque, which typically amounts to ~ 17 % (for a = 1.5, 
see Fig. [T2|) . This comparison indicates that the disk self-gravity does hardly change, if at 
all, the corotation torque. 



6. Concluding remarks 



The present work investigates the impact of the disk self-gravity on the type I migration. 
We show that the assumption customarily used in planet-disk calculations, namely a planet 
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freely migrating in a disk without self-gravity, can lead to a strong overestimate of the 
migration rate. We provide a simple evaluation of this overestimate (Fig. [5]). The drift rate 
can be overestimated by as much as a factor of two. Such a factor is inappropriate for the 
accurate calculation of migration rates, which is the main motivation of many recent studies 
of planet-disk interactions. The planet and the disk must therefore orbit within the same 
potential to yield unbiased estimates of the drift rate. Avoiding a spurious shift of resonances 
may be even more crucial in a non-barotropic situation. In this case, the corotation torqu e 



depends strongly upon the distance between orbit and corotation (IBaruteau fc Massetll2008l ). 
so that an ill-located corotation would yield meaningless drift rates. 

The inclusion of the disk self-gravity in our calculations confirms that the disk grav- 
ity accelerates t ype I migrati on. We solve the con tradiction between the statements of 
Nelson fc Benz and iPierens fc Hurg ( 120051 ) regarding the impact of the disk self- 



gravity on the migration rate. The increase of the differential Lindblad torque due to the 
disk gravity is typically one order of magnitude smaller than the spurious one induced by 
a planet freely migrating in a non self-gravitating disk. We provide a simple evaluation of 
this torque increase (Fig. [8]), which depends only on the Toomre parameter at the planet 
location, whatever the mass distribution of the whole disk. Furthermore, we argue that it 
can be entirely accounted for by a shift of the Lindblad resonances, and be modeled with 
an anisotropic pressure tensor. This model succeeds in reproducing the differential Lindblad 
torque of a self-gravitating calculation, but increases the corotation torque. This model 
enables us to conclude that there is no significant impact of the disk self-gravity on the 
corotation torque, in the linear regime. 

In a future work, we will extend our study beyond the linear regime. Preliminary 
calculations show that, regardless of the planet mass, the disk gravity speeds up migration. 
It would also be of interest to extend this study in three-dimensions. In the linear regime, 
we do not expect the torque relative increase due to the disk gravity to be altered in three- 
dimensions. However, three-dimensional calculations, involving the gas self-gravity, should 
be of considerable relevance for intermediate planet masses when a circumplanetary disk 
builds up, in particular to assess the frequency of type III migration. 

We thank the anonymous referee for a careful and insightful report. 



A. Expressions of Qr and g^p 



In this section, we give the expressions of the radial and azimuthal self-gravitating 
accelerations Qr and g^p, smoothed over the softening length e^g. We use the variables {u = 
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log(r/rinin), if), where Tmin denotes the inner edge radius of the grid. With this set of 
coordinates, gr{u,ip) reads 

gj.{u,(p) = — Ge""/^ / / Sr{u ,(p ) Kr{u — u ,(p — (fi ) du dip 



+ GJ:{u,ip)Kr{0,0)AuAip, (Al) 
where Sr and Kr are defined as 



n ^ N N «/2 1 / N 1 + - e "cos(v2) 

Sr{u, if) = S M, if) e I and K, U, Lp) = y— — / , T32 ul3/2 - 

|2(cosh(u) — cos(v9)) + B^e"}'^/^ 



(A2) 



In Eqs. ( lAip and (lA2p . G denotes the gravitational constant, Wmax = log(rmax/''"min) with 
'"max the outer edge radius of the grid, S is the disk surface density, Au and Ap are the mesh 
sizes, Kr{0,0) = 1/B and B = Ssg/r. Since egg oc r (see section 12.2. ip . B is uniform over 
the grid. The second term on the R.H.S. of Eq. flAip is an additional corrective term that 
ensures the absence of radial self-force. Similarly, g^{u, p) reads 

/•"max /'27r 

g^{u, p) = -Ge"^"/^ / / S^{u\ p) K^{u -u ,p) - p) du d(p\ (A3) 

Jo Jo 

with Sip and given by 

Spin, p>) = S(tx, p>) e^"/^ and K^(u, ^) = ''''^f^^_un2 ui3/2 - (^4) 

|2(cosh(u) — cos(v9)) + B^e^}-^/^ 

In the particular case where only the axisymmetric component of the disk self-gravity is 
accounted for, which involves the axisymmetric component of the disk surface density S(m) = 
(27r)^^ Jq'^ S(m, 'p)dp>, g^p cancels out and 



/* ^max 

griu) = -Ge-"/2 / Yriu) Kr{u - u) du + GT.{u)AuKr{^), (A5) 
Jo 



where Sr{u) = (27r) ^ J^^ Sr{u,(p)dp and Kr{u) = Jq"" Kr{u,(p)dp>. 



B. Relative difference of the torques between the free and fixed situations 

(without disk gravity) 

We denote by ST the difference of the one-sided Lindblad torques between the free and 
fixed cases. This difference can be written as 
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where x = r — Vp, 5x is the shift of the Lindblad resonances induce d by the free case, T is the 



m 



th 



Fourier component of the one-sided Lindblad torque (see e.g. I Ward! 119971 . or Eq. (|T4|) 



and xl is the location of the Lindblad resonances in the fixed situation: 



x, = -e 



2 ^T+e 



hr. 



(B2) 



with ^ = m/i, e = +1 for outer resonances, e = — 1 for inner resonances. Approximating the 
summation over m as an integral over ^, Eq. ( IBll) can be recast as 



= / {d^T/T) xTx6x d^. 



(B3) 

In Eq. (IB3p . T depends on x through the square of the forcing function \E', which i s usu ally 
approximated as a function of the Bessel functions Kq and Ki (see e.g. IWardI 119971 . or 
Eq. ( fTSl) ). Furthermore, "^( x.E) can be approximated as hrp/\x\C,, to within a numerical 
factor of the order of unity (lAbramowitz &: Stegunl Il972l ). Thus, T oc and d^T/T oc 
x^^. At the location of Lindblad resonances, given by Eq. (lB2p . this yields d^T/T oc eh^^. 

^. The shift Sx, which has same sign for inner and outer Lindblad 
5. The difference of the differential Lindblad t< 
obtained by summing Eq. (1B3P at inner and outer Lindblad resonances. 



Moreover, T oc eTjph 

resonances, scales with Sp. The difference of the differential Lindblad torques is eventually 



(5FiLR + 5FoLR oc T^ph ^ j (Tqlr - Tilr)^^ oc S^/i ^ x Sp/i ^ oc S^/i" 



(B4) 



Since the differential Lindblad torque scales with Ep/i ^, the relative difference of the dif- 
ferential Lindblad torques between the free and fixed cases scales with Sp/i~^, hence with 
{Qh)-\ 



C. Relative difference of the torques with and without disk gravity 

The calculation of the difference 5T of the one-sided Lindblad torques between the fully 
self-gravitating and non-gravitating situations is similar to the one derived in Appendix [Bl 
The difference 5T is given again by Eq. (IBSp , where 6x is this time the shift induced by the fsg 
situation. This shift has an opposite sign at inner and outer Lindblad resonances: 5x oc eSp 
(see the 5i?2 expression of PH05, or skip to Eq. (1121) where however x = (r — rp)/hrp). 
Furthermore, assumi ng that the expression of T given by Eq. (1141) can be applied for a 



self-gravitating disk (iGoldreich Sz Tremaind 119791 ). we still have dxT/T (x eh ^. Since the 



differential Lindblad torque scales with Ep/i ^, we find 

5FiLR + 5FoLR oc Sp/i-1 J (ToLR + Tilr)^^ oc X Sp/i-^ oc S^-^ (CI) 



-27- 



The relative difference of the differential Lindblad torques between the fsg and nog cases 
therefore scales with Ilph~^, hence with Q~^. 



D. Numerical and analytical shifts of Lindblad resonances induced by the disk 

self-gravity 

We studied in section [3.3.2l the relative difference of the torques between the fsg and nog 
situations. In particular, we find that our calculations with e = 0.1H{rp), which matches 
the mesh resolution in the planet vicinity, display a relative difference that is stronger than 
the one obtained with the estimate of PH05, which however does not involve a softening 
parameter. We give hereafter more insight into this result. 

We propose to evaluate for each azimuthal wavenumber m the shift of the Lindblad 
resonances induced by our fsg calculations, and compare it with its theoretical expression 
given by Eq. fll2p . This theoretical expression predicts that the shifts at inner and outer 
resonances are of opposite sign, their absolute value, which we denote by 6xth,m, being 
identical. We furthermore denote fenum.m the shift (in absolute value) inferred from our 
calculations, and (rfsg^m) the m^^ Fourier component of the inner (outer) Lindblad 

torque of a fsg calculation. We use similar notations for a nog calculation, and we drop 
hereafter the m subscripts for the sake of legibility. A first-order expansion yields 

-pi _ \ f) A-r anr\ — — f) r° A-r (^)^^ 

fsg nog ' '-'x'- nog "■^num cxiiu. ± jgg ± j^^g '-^x'- nog '-'•'^num- V / 

To estimate the quantities dxT\^og ^i^F^^g, we performed an additional nog calculation, 
mentioned as nogo calculation, for which we imposed a slight, known shift of the resonances. 
This was done by fixing the planet's angular velocity at Vtp — 5VLp, with 5VLp = 10~^ flp. This 
slight decrease of the planet's angular velocity, with respect to the nog situation, implies an 
outward shift of inner and outer Lindblad resonances that reads 6x0 = {26Qp) / {3hQp) , ex- 
pression that is independent of m. With similar notations as before for the nogo calculation, 
and using again a first-order expansion, we have 

r* — r* -L f) nnH r° — r° -l ft r° (r)9\ 

^ nogo ~ nog ' nog "-^o aifj. ± ^^g^ — ± T nog '^-''O- \^^) 

Combining Eqs. (IDip and (]D2p . we are finally left with 

SXnu^ = 1'^' J''\ ^ X SXo. (D3) 

V nogo ' nogo/ V nog ' nog/ 



We plot in Fig. [TUthe ratio Sx^^^^/Sxth as a function of the azimuthal wavenumber m. 



for Sp = 2 X 10 ^ (Q ~ 8). We first comment that the ratio is negative for m < 6, positive 
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beyond, with a divergent behavior at the transition. We checked that this behavior is caused 
by a change of sign of the denominatoi0 of Eq. flD3p . which is negative for m < 6 and 
positive beyond. Furthermore, the ratio Sxnum/ ^Xt^ is significantly greater than unity for m 
ranging from ~ 7 to ~ 20, that is for the dominant Lindblad resonances. Differently stated, 
the dominant Lindblad resonances are more shifted by our calculations than analytically 
expected by PH05, which explains why the torque enhancement is more important with our 
calculations. 

Beyond, the ratio is close to unity for a rather large range of high m- values. This confirms 
that for high values of m the WKB approximation yields analytical estimates that are in 
good agreement with the results of numerical simulations. However, since our calculations 
involve a softening parameter, the ratio does not converge when increasing m, and slowly 
tends to zero. We checked that the value of m for which the ratio becomes lower than unity 
increases when decreasing the softening length. This explains why the torque enhancement 
is increasingly important at smaller softening length, as inferred from Fig. [TOl 
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Fig. 2. — Left: Relative difference of the torques obtained without self-gravity (Fnog), and 
with self-gravity (Ffsg), for different grid resolutions (see text). Right: Torque exerted on a 
Mp — 5 X 10~^ planet mass, and on the system {gas-|-planet}. Torques are depicted for a 
calculation without self-gravity {long-dashed and dash-dotted curves), and with self-gravity 
{solid and dotted curves). 
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Fig. 3. — Specific torque 7 exerted on a Mp = 5 x 10~^M* planet mass by a non self- 
gravitating disk, with a h = 5 % aspect ratio. Diamonds refer to the fixed case (the planet is 
held on a fixed circular orbit, with a strictly Keplerian angular velocity) while asterisks refer 
to the free case (the planet freely evolves in the disk, the planet's angular velocity is greater 
than Keplerian). The solid line corresponds to a proportional fit of the fixed case data, 
and shows the excellent agreement between our results of calculations and the expectation 
7 oc Sp in the fixed case. The two dott ed lines depict the two- and three-dimensional 
analytical estimates of iTanaka et al.l (120021 ). 
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Fig. 4. — Location of two Lindblad resonances in the fixed case (left panel) and in the free case 
(right panel): the ILR of m = 6 (^^ilr = 6/5 ilp), and the OLR of m = 5 (f2oLR = 5/6 Qp). 
The disk's rotation profile f2(r) is depicted without self-gravity (sohd curve) and with self- 
gravity (dashed curve, right panel). In the latter case, gr{r) is given by a calculation with 
T,p — 5 X 10~^, a value exaggerated to improve legibility. Note also that the pressure 
buffer has been discarded in both profiles, for the sake of simphcity. The vertical arrow 
at r = 1 indicates the planet location, it reaches the upper curve in the free case (right 
panel) since the planet feels the disk gravity. The ILR and OLR are located, respectively, at 
Hlr = ^^^(^ilr) and tolr = ^~^(^olr)- The nominal position of the resonances (that of 
the fixed case) is indicated by light gray dash-dotted lines on the right panel to appreciate 
their shift, highlighted by a horizontal arrow. 
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Fig. 5. — Relative difference of tlie torques between tlie free and fixed situations, as a function 
of Qh [see text and Eq. ([7])]. 
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Fig. 6. — Specific torque on a Mp = 5 x lO^^M* planet mass, obtained with axisymmetric 
and fully self-gravitating calculations, with a h = 5 % disk aspect ratio. Torques achieved 
without disk gravity (see section 13. ip are also displayed, for comparison. 
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Fig. 7. — Same as Fig. HI except that we examine the shift of the LR when the disk is self- 
gravitating (its rotation profile is now the solid, upper curve). In the left panel, the planet 
does not feel the disk gravity: the frequency of the planet, and therefore that of the LR, is 
the same as in Fig. H^. In the right panel, both the planet and the disk feel the disk gravity: 
the frequencies of the planet and of the LR are those of Fig. H^. 
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Fig. 8. — Relative difference of the torques obtained with the fully self-gravitating calcu- 
lations (7fsg) and the calculations without disk gravity (7nog), as a function of the Toomre 
parameter Q at the planet location. 
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Fig. 9. — Specific torque variation with time, with and without disk gravity. In each case, 
two situations are depicted: the fixed case (the planet is on a fixed orbit with the appropriate 
angular velocity, see Table [2]) and the free case (the planet is free to migrate in the disk). 
Except the self-gravitating calculation with a free planet, the calculations are those of Figs. [3] 
and E] for Ep = 2 x 10^^. When the planet in on a free orbit without self-gravity, the torque 
oscillates with a large amplitude. This is due to the slight increase of the planet frequency: 
^p(^p); which is initially strictly Keplerian, is given by Eq. ([5]) during its time evolution. 
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Fig. 10. — Relative difference of the torques between with the fsg situation (7fsg) and the 
nog situation (7nog) as a function of the disk surface density Ep. We compare the results of 
our calculations (each symbol refers to a different value of the softening length e) with the 
analytical results of PH05 (dashed curve), and our analytical estimate (dash-dotted curve, 
see text and Eq. (ITTl) ). The vertical dotted lines display different values of the Toomre 
parameter at the planet location. 
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Fig. 11. — Specific torque exerted on a Mp — 10~^M* planet mass, as a function of the 
Toomre parameter Q at the planet location. We display the torques obtained with asg 
calculations (plus signs), fsg calculations (asterisks) and asg+ap calculations (diamonds). 
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Fig. 12. — Specific torque as a function of the planet to primary mass ratio. Calculations 
obtained with the anisotropic pressure model (asg+ap) are compared with axisymmetric 
self-gravitating calculations (asg) and fully self-gravitating calculations (fsg). The close-up 
displays the relative difference of the torques between the fsg and asg+ap situations. For all 
these calculations, Q — S at the planet location (the disk mass is ~ 0.024M*). 
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Fig. 13. — Specific torque as a function of tlie planet to primary mass ratio, for a flat initial 
surface density profile. The square corresponds to an additional nog calculation performed 
with g = 5 X 10~^. The close-up displays the relative difference of the corotation torques 
between the asg and fsg situations (see text and Eq. ( 126|) ). 
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Fig. 14. — Ratio of 5xnunu the shift of Lindblad resonances obtained with our fsg calculations 
(see Eq. flD3l) ). and of 5xth, its analytically expected value (see Eq. f|T2l) ). 



